In a simplest version of radiation dose optimization, a dose matrix A εn×m indicates a radiation target dose delivered by a set of m beams in n voxels. Here, the voxels represent small cubes (mm3) in a volume of tissue. A vector of nonnegative beam weights xεm linearly combines the beams to yield a vectory=Axεn of cumulative voxel radiation dosages. Typically, the dose matrix A is relatively tall (n>>m), and sparse. The problem is to determine beam weights that guarantee no underdoses and minimal overdoses, given a target dosage pattern bεn. This is a linear program (LP):
                                                        min              x                        ⁢                                                                                                                        A                      ⁢                                                                                          ⁢                      x                                        -                    b                                                                    1                            ⁢                                                          ⁢              such              ⁢                                                          ⁢              that              ⁢                                                          ⁢              A              ⁢                                                          ⁢              x                                ≥          b                ,                  x          ≥          0.                                    (        1.1        )            
A small problem (m˜103 beams and n˜2m voxels) can currently be solved in roughly one second on a 2 GHz Pentium 4 processor. However, the LP solution time generally scales as the cube of the problem size, and so that method is presently impractical for practical clinical problems with much larger sizes.
LP and second-order cone (SOCP) methods are attractive because they lend themselves to formulations that minimize dose delivery errors. However, problems are limited to just a few thousand variables and constraints, and require hours to solve. For that reason, the problem prefers an L2 norm in an objective function, i.e., minimize the sum squared error while disallowing underdoses. This is a least-squares with inequalities problem (LSI):
                                                        min              x                        ⁢                                                                                                                        A                      ⁢                                                                                          ⁢                      x                                        -                    b                                                                    2                2                            ⁢                                                          ⁢              such              ⁢                                                          ⁢              that              ⁢                                                          ⁢              A              ⁢                                                          ⁢              x                                ≥          b                ,                  x          ≥          0.                                    (        1.2        )            
The LSI can be justified as optimizing an upper bound on the LP. Minimizing sum-squared overdose gives solutions that are slightly smoother, but not as safe because the overdose is larger. This LSI is a particularly well-behaved instance of a quadratic program (QP), assuming the matrix A has a full row rank, and the QP is strictly convex. Herein, LSI and QP are used interchangeably.
Fast Approximate Methods
In the prior art, it is common to use methods that avoid the underdosing constraint and simply seek to reduce a sum squared error∥Ax−b∥22.or a weighted sum squared error∥W(Ax−b)∥22.
Due to the size of the problem, solutions are computed iteratively via gradient, conjugate gradient, or quasi-Newton methods, or preferably, a limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method. To avoid physically impossible solutions, those methods are usually modified so any beam weight xi that is negative is reset to zero. That is not guaranteed to find an optimum or avoid underdoses, but appears to work well enough for clinical use.
More principled methods, such as gradient projection, can guarantee non-negativity and optimality, but increase processing time and are thus slower. For this reason, parallel solutions on graphic processor units (GPU) are sometimes used that solve small problems in a matter of seconds.